Creating codebooks
The authors provide two datasets, but no further documentation -
thus, we started by attempting to automatically generate some codebooks
to understand the data
library(haven)
library(gt)
library(patchwork)
library(modelsummary)
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
## backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
##
## Revert to `kableExtra` for one session:
##
## options(modelsummary_factory_default = 'kableExtra')
## options(modelsummary_factory_latex = 'kableExtra')
## options(modelsummary_factory_html = 'kableExtra')
##
## Silence this message forever:
##
## config_modelsummary(startup_message = FALSE)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(broom)
library(metafor)
## Loading required package: metadat
## Loading required package: numDeriv
##
## Loading the 'metafor' package (version 4.6-0). For an
## introduction to the package please type: help(metafor)
library(broom.mixed)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
dat1 <- read_dta("reproduction_materials/combined data OSF.dta")
Table 1 (excluding Study 6)
Likely typo in article - Table 1 states N Study 1 as 101, but data is
102, as is in note to Fig 2 (noone was excluded in that study, so
changing sample size is odd) Missing condition information for 1 out of
102 - so should have been excluded
dat1 %>% count(study, sample) %>% group_by(study) %>%
summarise(N_total = sum(n), N_analytic = n[sample == 1]) %>% left_join(
dat1 %>% group_by(study) %>%
summarise(
Real_effort = any(!is.na(numbersearches)),
Stated_effort = any(!is.na(tokens)),
Bonus_expectations = (any(!is.na(s2_tokens_0
)) | any(!is.na(muchlower_score))))
) %>% gt()
## Joining with `by = join_by(study)`
| study |
N_total |
N_analytic |
Real_effort |
Stated_effort |
Bonus_expectations |
| 1 |
102 |
102 |
FALSE |
TRUE |
FALSE |
| 2 |
103 |
96 |
FALSE |
TRUE |
TRUE |
| 3 |
68 |
68 |
TRUE |
FALSE |
TRUE |
| 4 |
306 |
258 |
TRUE |
FALSE |
TRUE |
| 5 |
605 |
518 |
TRUE |
FALSE |
TRUE |
# Exclude excluded participants
dat1 <- dat1 %>% filter(sample == 1)
Figure 2
s2 <- dat1 %>%
filter(study == 2) %>%
select(starts_with("s2_tokens_"), cond) %>%
pivot_longer(
cols = starts_with("s2_tokens"),
names_to = "tokens",
values_to = "expectation",
names_prefix = "s2_tokens_"
) %>% filter(tokens != "mean") %>%
mutate(tokens = as.numeric(tokens)) %>%
group_by(tokens, cond) %>%
summarise(mean = mean(expectation, na.rm = TRUE),
se = sd(expectation, na.rm = TRUE) / sqrt(sum(!is.na(expectation))),
ci_lower = mean - 1.96 * se,
ci_upper = mean + 1.96 * se) %>%
ggplot(aes(x = tokens, y = mean, color = haven::as_factor(cond))) +
geom_line() +
scale_y_continuous(limits = c(1, 5)) +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
labs(x = "Tokens contributed", y = "Subjective likelihood of bonus", col = "") +
scale_color_manual(values = c("Control" = "blue", "Advantage" = "green", "Disadvantage" = "orange")) +
theme_minimal()
## `summarise()` has grouped output by 'tokens'. You can override using the
## `.groups` argument.
s3_5 <- dat1 %>%
filter(study > 2) %>%
select(matches("_score"), cond) %>%
pivot_longer(
cols = matches("_score"),
names_to = c("performance", NA),
values_to = "expectation",
names_sep = "_"
) %>%
mutate(performance = fct_recode(performance,
"Much Lower" = "muchlower",
"Lower" = "lower",
"About Same" = "same",
"Higher" = "higher",
"Much Higher" = "muchhigher") %>%
factor(levels = c("Much Lower", "Lower", "About Same", "Higher", "Much Higher"))) %>%
group_by(performance, cond) %>%
summarise(mean = mean(expectation, na.rm = TRUE),
se = sd(expectation, na.rm = TRUE) / sqrt(sum(!is.na(expectation))),
ci_lower = mean - 1.96 * se,
ci_upper = mean + 1.96 * se) %>%
ggplot(aes(x = performance, y = mean, color = haven::as_factor(cond), group = haven::as_factor(cond), )) +
geom_line() +
scale_y_continuous(limits = c(0, 100)) +
geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), width = 0.2) +
labs(x = "Performance relative to other employee(s)", y = "Subjective likelihood of bonus (%)", col = "") +
scale_color_manual(values = c("Control" = "blue", "Advantage" = "green", "Disadvantage" = "orange")) +
theme_minimal()
## `summarise()` has grouped output by 'performance'. You can override using the
## `.groups` argument.
p1 <- s2 / s3_5
p1
Control group mean misreported in text, is mean of disadvantaged
group
Means
s2_means <- dat1 %>%
filter(study == 2) %>%
rowid_to_column() %>%
select(-tokens) %>%
pivot_longer(
cols = starts_with("s2_tokens"),
names_to = "tokens",
values_to = "expectation",
names_prefix = "s2_tokens_"
) %>%
group_by(rowid) %>%
summarise(cond = as_factor(cond)[1],
m_effort = mean(as.numeric(tokens), na.rm = TRUE), # Should be fixed, just code check
m_expectation = mean(expectation, na.rm = TRUE),
sd_expectation = sd(expectation, na.rm = TRUE)
)
## Warning: There were 96 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `m_effort = mean(as.numeric(tokens), na.rm = TRUE)`.
## ℹ In group 1: `rowid = 1`.
## Caused by warning in `mean()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 95 remaining warnings.
mod_s2_avg <- lm(m_expectation ~ cond, s2_means)
modelsummary(mod_s2_avg, estimate = "{estimate} [{conf.low}, {conf.high}]",
statistic = "({std.error}), p = {p.value}",
gof_map = NA,
notes = as.character(glue::glue("df = {modelsummary::glance(mod_s2_avg)$df.residual}")),
title = "Average expectation in Study 2")
tinytable_l28mme817efhw05euxa1
Average expectation in Study 2
| |
(1) |
| df = 93 |
| (Intercept) |
2.961 [2.776, 3.146] |
| |
(0.093), p = <0.001 |
| condAdvantage |
0.426 [0.172, 0.680] |
| |
(0.128), p = 0.001 |
| condDisadvantage |
-0.586 [-0.844, -0.328] |
| |
(0.130), p = <0.001 |
s3_5_means <- dat1 %>%
filter(study > 2) %>%
rowid_to_column() %>%
select(-tokens) %>%
pivot_longer(
cols = matches("_score"),
names_to = c("performance", NA),
values_to = "expectation",
names_sep = "_"
) %>%
group_by(rowid) %>%
summarise(study = study[1],
cond = as_factor(cond)[1],
effort_ratings = n(), # Should be fixed, just code check
m_expectation = mean(expectation, na.rm = TRUE),
sd_expectation = sd(expectation, na.rm = TRUE)
)
mod_s3_5_avg <- lmer(m_expectation ~ cond + (1 | study), s3_5_means)
modelsummary(mod_s3_5_avg,
estimate = "{estimate} [{conf.low}, {conf.high}]",
statistic = "({std.error}), p = {p.value}",
title = "Average expectation in pooled studies 3-5")
tinytable_7rp384wkw4uuowymdy01
Average expectation in pooled studies 3-5
| |
(1) |
| (Intercept) |
49.363 [44.396, 54.330] |
| |
(2.531), p = <0.001 |
| condAdvantage |
15.203 [13.007, 17.399] |
| |
(1.119), p = <0.001 |
| condDisadvantage |
-14.225 [-16.413, -12.037] |
| |
(1.115), p = <0.001 |
| SD (Intercept study) |
4.107 |
| |
, p = |
| SD (Observations) |
13.301 |
| |
, p = |
| Num.Obs. |
844 |
| R2 Marg. |
0.423 |
| R2 Cond. |
0.473 |
| AIC |
6772.1 |
| BIC |
6795.8 |
| ICC |
0.1 |
| RMSE |
13.26 |
Slopes
s2_expectations <- dat1 %>%
filter(study == 2) %>%
rowid_to_column() %>%
select(-tokens) %>%
pivot_longer(
cols = starts_with("s2_tokens"),
names_to = "tokens",
values_to = "expectation",
names_prefix = "s2_tokens_"
) %>% filter(tokens != "mean") %>%
mutate(tokens = as.numeric(tokens), cond = as_factor(cond))
s2_expectations <- s2_expectations %>% group_by(cond, rowid) %>%
mutate(delta = expectation - lag(expectation)) %>%
summarise(delta = mean(delta, na.rm = TRUE))
## `summarise()` has grouped output by 'cond'. You can override using the
## `.groups` argument.
mod_s2_slope <- lm(delta ~ cond, s2_expectations) %>% summary()
modelsummary(mod_s2_slope, estimate = "{estimate} [{conf.low}, {conf.high}]",
statistic = "({std.error}), p = {p.value}",
gof_map = NA,
notes = as.character(glue::glue("df = {modelsummary::glance(mod_s2_slope)$df.residual}")),
title = "Expected reward sensitivity in Study 2")
tinytable_h1yimrxa5jam6ka81r5o
Expected reward sensitivity in Study 2
| |
(1) |
| df = 93 |
| (Intercept) |
0.753 [0.660, 0.847] |
| |
(0.047), p = <0.001 |
| condAdvantage |
-0.165 [-0.294, -0.037] |
| |
(0.065), p = 0.012 |
| condDisadvantage |
-0.253 [-0.384, -0.123] |
| |
(0.066), p = <0.001 |
s3_5_expectations <- dat1 %>%
filter(study > 2) %>%
rowid_to_column() %>%
select(-tokens) %>%
pivot_longer(
cols = matches("_score"),
names_to = c("performance", NA),
values_to = "expectation",
names_sep = "_"
) %>% mutate(performance = fct_recode(performance,
"Much Lower" = "muchlower",
"Lower" = "lower",
"About Same" = "same",
"Higher" = "higher",
"Much Higher" = "muchhigher") %>%
factor(levels = c("Much Lower", "Lower", "About Same", "Higher", "Much Higher")),
cond = as_factor(cond)) %>%
group_by(study, cond, rowid) %>%
mutate(delta = expectation - lag(expectation)) %>%
summarise(delta = mean(delta, na.rm = TRUE))
## `summarise()` has grouped output by 'study', 'cond'. You can override using the
## `.groups` argument.
mod_s3_5_delta <- lmer(delta ~ cond + (1 | study), s3_5_expectations)
## boundary (singular) fit: see help('isSingular')
modelsummary(mod_s3_5_delta,
estimate = "{estimate} [{conf.low}, {conf.high}]",
statistic = "({std.error}), p = {p.value}",
title = "Expected reward sensitivity in pooled studies 3-5")
tinytable_sco7nb3uh7h213j81zk7
Expected reward sensitivity in pooled studies 3-5
| |
(1) |
| (Intercept) |
22.178 [21.426, 22.930] |
| |
(0.383), p = <0.001 |
| condAdvantage |
-6.124 [-7.205, -5.044] |
| |
(0.551), p = <0.001 |
| condDisadvantage |
-7.995 [-9.073, -6.918] |
| |
(0.549), p = <0.001 |
| SD (Intercept study) |
0.000 |
| |
, p = |
| SD (Observations) |
6.545 |
| |
, p = |
| Num.Obs. |
843 |
| R2 Marg. |
0.217 |
| AIC |
5567.0 |
| BIC |
5590.7 |
| RMSE |
6.53 |
Figure 2 (& table)
# Calculate standardised effects
summaries <- dat1 %>%
filter(!is.na(cond)) %>%
group_by(study) %>%
mutate(z_effort = (coalesce(log(numbersearches), tokens) - mean(coalesce(log(numbersearches), tokens), na.rm = TRUE))/sd(coalesce(log(numbersearches), tokens), na.rm = TRUE)) %>%
group_by(study, cond) %>%
summarise(cond = as_factor(cond)[1],
m_effort = mean(z_effort, na.rm = TRUE),
sd_effort = sd(z_effort, na.rm = TRUE),
N = n(),
se = sd_effort / sqrt(N),
.groups = "drop") %>%
mutate(conf.low = m_effort - 1.96 * se,
conf.high = m_effort + 1.96 * se,
vi = se^2)
# Add p-values per study
p_values <- dat1 %>%
filter(!is.na(cond)) %>%
group_by(study) %>%
mutate(z_effort = (coalesce(log(numbersearches), tokens) -
mean(coalesce(log(numbersearches), tokens), na.rm = TRUE)) /
sd(coalesce(log(numbersearches), tokens), na.rm = TRUE),
cond = cond %>% as_factor() %>% relevel(ref = "Control")) %>%
group_by(study) %>%
do({
data <- .
model <- lm(z_effort ~ cond, data = data)
tidy(model) %>%
filter(term != "(Intercept)") %>% # Remove the intercept term
mutate(study = unique(data$study))
}) %>%
transmute(study, cond = term %>% str_remove("cond"), p_value = p.value) %>%
ungroup()
summaries <- summaries %>% left_join(p_values)
## Joining with `by = join_by(study, cond)`
# Run meta-analyses
## Model without an intercept to get estimates and standard errors for each condition
model_no_intercept <- rma.mv(yi = m_effort, V = vi, mods = ~ factor(cond) - 1, random = ~ 1 | study, data = summaries)
estimates_no_intercept <- coef(summary(model_no_intercept))[, "estimate"]
se_no_intercept <- coef(summary(model_no_intercept))[, "se"]
## Model with condition as a moderator to get p-values comparing each condition to control
summaries$cond <- relevel(factor(summaries$cond), ref = "Control")
model_with_intercept <- rma.mv(yi = m_effort, V = vi, mods = ~ cond, random = ~ 1 | study, data = summaries)
summary_with_intercept <- summary(model_with_intercept)
## Extract p-values from the model with intercept
p_values <- coef(summary_with_intercept)[-1, "pval"]
p_values <- c(NA, p_values) # Control comparison is NA
meta_results <- tibble(
study = "meta",
cond = levels(factor(summaries$cond)),
m_effort = estimates_no_intercept,
se = se_no_intercept,
p_value = p_values,
conf.low = m_effort - 1.96 * se,
conf.high = m_effort + 1.96 * se
)
fig2_res <- meta_results %>% bind_rows(summaries %>% mutate(study = as.character(study)), .)
# Function to calculate t-test p-value using summary statistics
fig2_res %>% select(-vi) %>%
gt() %>%
fmt_number(columns = c(m_effort, conf.low, conf.high, se, sd_effort), decimals = 3) %>%
fmt(
columns = p_value,
fns = \(p) timesaveR::fmt_p(p, equal_sign = FALSE)
)
| study |
cond |
m_effort |
sd_effort |
N |
se |
conf.low |
conf.high |
p_value |
| 1 |
Control |
0.437 |
0.971 |
34 |
0.167 |
0.111 |
0.764 |
NA |
| 1 |
Advantage |
0.026 |
0.982 |
32 |
0.174 |
−0.314 |
0.366 |
.079 |
| 1 |
Disadvantage |
−0.448 |
0.866 |
35 |
0.146 |
−0.735 |
−0.162 |
< .001 |
| 2 |
Control |
0.365 |
0.873 |
30 |
0.159 |
0.053 |
0.678 |
NA |
| 2 |
Advantage |
0.092 |
0.933 |
34 |
0.160 |
−0.222 |
0.406 |
.256 |
| 2 |
Disadvantage |
−0.440 |
1.042 |
32 |
0.184 |
−0.801 |
−0.079 |
.001 |
| 3 |
Control |
0.174 |
1.073 |
23 |
0.224 |
−0.265 |
0.612 |
NA |
| 3 |
Advantage |
0.013 |
1.008 |
22 |
0.215 |
−0.408 |
0.434 |
.593 |
| 3 |
Disadvantage |
−0.186 |
0.925 |
23 |
0.193 |
−0.564 |
0.192 |
.228 |
| 4 |
Control |
0.271 |
0.842 |
91 |
0.088 |
0.098 |
0.444 |
NA |
| 4 |
Advantage |
−0.025 |
0.942 |
85 |
0.102 |
−0.226 |
0.175 |
.046 |
| 4 |
Disadvantage |
−0.274 |
1.143 |
82 |
0.126 |
−0.522 |
−0.027 |
< .001 |
| 5 |
Control |
0.152 |
0.895 |
178 |
0.067 |
0.020 |
0.283 |
NA |
| 5 |
Advantage |
0.029 |
0.980 |
167 |
0.076 |
−0.120 |
0.177 |
.250 |
| 5 |
Disadvantage |
−0.184 |
1.093 |
173 |
0.083 |
−0.347 |
−0.021 |
.002 |
| meta |
Control |
0.020 |
NA |
NA |
0.052 |
−0.083 |
0.123 |
NA |
| meta |
Advantage |
0.229 |
NA |
NA |
0.047 |
0.136 |
0.322 |
.003 |
| meta |
Disadvantage |
−0.266 |
NA |
NA |
0.057 |
−0.378 |
−0.155 |
< .001 |
# Convert p-values to a readable format for annotations
fig2_res <- fig2_res %>%
mutate(p_label = ifelse(is.na(p_value), "", paste0("p ", timesaveR::fmt_p(p_value, equal_sign = TRUE, digits = 3))))
# Create the plot using ggplot2
# Reorder the conditions
fig2_res$cond <- factor(fig2_res$cond, levels = c("Control", "Advantage", "Disadvantage"))
# Create the plot using ggplot2
ggplot(fig2_res, aes(x = as.factor(study), y = m_effort, color = cond, group = cond)) +
geom_point(position = position_dodge(width = 0.7), size = 3, shape = 16) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(width = 0.7)) +
geom_text(aes(label = p_label, y = conf.low), position = position_dodge(width = 0.7), vjust = 1.4, size = 3, color = "black") +
labs(x = "Study", y = "Standardized mean", color = NULL) +
theme_minimal() +
theme(legend.position = "bottom", legend.title = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_color_manual(values = c("Control" = "blue", "Advantage" = "darkgreen", "Disadvantage" = "darkorange")) +
scale_x_discrete(labels = c("1" = "Study 1", "2" = "Study 2", "3" = "Study 3", "4" = "Study 4",
"5" = "Study 5", "meta" = "Mini meta-analysis")) +
guides(color = guide_legend(override.aes = list(shape = 16))) +
expand_limits(y = c(min(fig2_res$conf.low) - 0.1, max(fig2_res$conf.high)))
# Mediation
# Study 2
med_data <- dat1 %>%
filter(!is.na(cond)) %>%
group_by(study) %>%
mutate(z_perf_slope = coalesce(scale(s2_performance_slope) %>% as.numeric(),
scale(performance_slope) %>% as.numeric()),
z_numbersearches = scale(numbersearches) %>% as.numeric(),
z_tokens = scale(tokens) %>% as.numeric(),
cond = as_factor(cond)) %>%
ungroup() %>%
filter(!is.na(z_perf_slope)) %>%
mutate(timesaveR::dummy_code(cond))
Regression results for Study 2
slope_mod_s2 <- med_data %>% filter(study == 2) %>% lm(z_tokens ~ z_perf_slope, .)
slope_mod_s2 %>% tidy(conf.int = TRUE) %>%
mutate(df = slope_mod_s2$df.residual) %>%
filter(term != "(Intercept)") %>%
glimpse()
## Rows: 1
## Columns: 8
## $ term <chr> "z_perf_slope"
## $ estimate <dbl> 0.3884817
## $ std.error <dbl> 0.09504098
## $ statistic <dbl> 4.087518
## $ p.value <dbl> 9.181902e-05
## $ conf.low <dbl> 0.1997756
## $ conf.high <dbl> 0.5771878
## $ df <int> 94
Regression results for Study 3-5
library(lmerTest)
slope_mod_s3_5 <- med_data %>% filter(study > 2) %>%
lmer(z_numbersearches ~ z_perf_slope + (1 | study), .)
## boundary (singular) fit: see help('isSingular')
slope_mod_s3_5 %>% summary()
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_numbersearches ~ z_perf_slope + (1 | study)
## Data: .
##
## REML criterion at convergence: 2375.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6412 -0.7776 -0.2550 0.6303 2.7683
##
## Random effects:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.0000 0.0000
## Residual 0.9714 0.9856
## Number of obs: 843, groups: study, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.167e-03 3.395e-02 8.410e+02 0.034 0.973
## z_perf_slope 1.659e-01 3.401e-02 8.410e+02 4.877 1.29e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## z_perf_slop 0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
Mediation results
Note that N in Fig 3 is reported as 844, so that it corresponds to
these results - even though it is not stated that the Fig only contains
3-5 rather than 2-5. Pattern for S2 looks very similar, though effects
are larger (in line with regression results).
library(lavaan)
## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.
med_data <- med_data %>% mutate(z_effort = coalesce(z_tokens, z_numbersearches))
mod_full <- ("
# Path from condition (cond) to mediator (z_perf_slope)
z_perf_slope ~ a1 * cond_advantage + a2 * cond_disadvantage
# Path from condition (cond) to outcome (z_effort)
z_effort ~ Advantage * cond_advantage + Disadvantage * cond_disadvantage + b * z_perf_slope
")
mod_reduced <- ("
z_effort ~ Advantage * cond_advantage + Disadvantage * cond_disadvantage
")
med_res <- sem(mod_full, data = med_data %>% filter(study > 2), cluster = "study") %>%
tidy(conf.int = TRUE) %>%
filter(str_detect(label, "Advantage|Disadvantage")) %>%
transmute(model = "full", cond = label, study = "3-5", estimate, conf.low, conf.high, p.value) %>%
bind_rows(
sem(mod_reduced, data = med_data %>% filter(study > 2), cluster = "study") %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "reduced", cond = label, study = "3-5", estimate, conf.low, conf.high, p.value)
) %>%
bind_rows(
sem(mod_full, data = med_data %>% filter(study == 2)) %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "full", cond = label, study = "2", estimate, conf.low, conf.high, p.value) %>% bind_rows(
sem(mod_reduced, data = med_data %>% filter(study == 2)) %>% tidy(conf.int = TRUE) %>% filter(str_detect(label, "Advantage|Disadvantage")) %>% transmute(model = "reduced", cond = label, study = "2", estimate, conf.low, conf.high, p.value)
)
) %>%
mutate(p_label = paste0("p ", timesaveR::fmt_p(p.value, equal_sign = TRUE, digits = 3)))
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
## The variance-covariance matrix of the estimated parameters (vcov)
## does not appear to be positive definite! The smallest eigenvalue
## (= -2.036678e-18) is smaller than zero. This may be a symptom that
## the model is not identified.
## Warning in lav_model_vcov(lavmodel = lavmodel, lavsamplestats = lavsamplestats, : lavaan WARNING:
## The variance-covariance matrix of the estimated parameters (vcov)
## does not appear to be positive definite! The smallest eigenvalue
## (= 6.488542e-20) is close to zero. This may be a symptom that the
## model is not identified.
med_res %>%
gt() %>%
fmt_number(columns = c(estimate, conf.low, conf.high), decimals = 3) %>%
fmt_number(columns = p.value, decimals = 4)
| model |
cond |
study |
estimate |
conf.low |
conf.high |
p.value |
p_label |
| full |
Advantage |
3-5 |
−0.035 |
−0.099 |
0.029 |
0.2859 |
p = .286 |
| full |
Disadvantage |
3-5 |
−0.154 |
−0.215 |
−0.092 |
0.0000 |
p < .001 |
| reduced |
Advantage |
3-5 |
−0.151 |
−0.306 |
0.004 |
0.0563 |
p = .056 |
| reduced |
Disadvantage |
3-5 |
−0.306 |
−0.458 |
−0.153 |
0.0001 |
p < .001 |
| full |
Advantage |
2 |
−0.083 |
−0.536 |
0.369 |
0.7181 |
p = .718 |
| full |
Disadvantage |
2 |
−0.514 |
−0.993 |
−0.036 |
0.0351 |
p = .035 |
| reduced |
Advantage |
2 |
−0.273 |
−0.734 |
0.188 |
0.2455 |
p = .245 |
| reduced |
Disadvantage |
2 |
−0.805 |
−1.273 |
−0.338 |
0.0007 |
p < .001 |
med_res %>%
group_by(study, cond) %>%
summarise(
share_mediated = 1 - (estimate[1] / estimate[2]), order = paste(model[1], model[2]),
.groups = "drop_last"
) %>%
dplyr::select(-order) %>%
gt() %>%
fmt_percent(columns = share_mediated, decimals = 1)
| cond |
share_mediated |
| 2 |
| Advantage |
69.4% |
| Disadvantage |
36.1% |
| 3-5 |
| Advantage |
76.8% |
| Disadvantage |
49.7% |
med_res %>%
mutate(model = factor(model, levels = c("reduced", "full"))) %>%
ggplot(aes(x = cond, y = estimate, group = model, shape = model, color = model)) +
geom_point(size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
width = 0.2, position = position_dodge(width = 0.4)
) +
geom_text(aes(y = conf.low - 0.05, label = p_label),
color = "black",
position = position_dodge(width = 0.4),
vjust = 1, size = 3
) +
facet_wrap(~study, ncol = 2, labeller = labeller(study = c("2" = "Study 2", "3-5" = "Study 3-5"))) +
scale_shape_manual(values = c(16, 1), labels = c("Reduced model", "Full model")) +
scale_color_manual(values = c("black", "black"), labels = c("Reduced model", "Full model")) +
labs(x = "", y = "Standardized effect", shape = "", color = "") +
scale_y_continuous(limits = c(-1.4, 0.5), breaks = round(seq(-1.4, 0.4, by = 0.2), 2)) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.title.y = element_text(margin = margin(r = 10)),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank()
) +
geom_hline(yintercept = 0, linetype = "solid", color = "black")
